Introduction

The response of ice caps on Severnaya Zemlya to climate change is poorly constrained due to their remote location within the Russian High Arctic (Dowdeswell et al., 2002, in press). Bassford et al. (2006) present a model of the mass balance of the Vavilov Ice Cap, October Revolution Island, Severnaya Zemlya (Fig. 1). The model successfully simulates the measured mass balance and evolution of englacial temperature. The results show that the mass balance is heavily dependent on processes leading to the refreezing of melted snow and ice. While surface processes are critical to the overall ice cap behavior, so too is the flow of ice.

Observations of the geometry of the Vavilov Ice Cap between 1952 and 1985 (Figs. 1 and 2a and b) allow us to hypothesize that the ice cap is not in equilibrium with the present climate. Barkov et al. (1992) compared maps of the Vavilov Ice Cap constructed from aerial photographs acquired in 1952 and 1985, in order to assess changes in the configuration of the ice cap. Cartographic analysis indicates that during this period the southern margins of the ice cap advanced by 150–450 m, while the northern and eastern margins remained stable or retreated slightly (Fig. 3). In contrast to a general retreat of the western margins, a 9 km section of the ice cap bordering the coast advanced by up to 400 m. Overall, an area of 11.1 km² was deglaciated between 1952 and 1985, compensated by a newly glaciated area of 14.6 km², with a resulting net increase of 3.5 km² in the area of the ice cap (Barkov et al., 1992). These observations are supported by recent ground-based measurements which indicate that, during the period 1984–1996, the northwestern margin of the Vavilov Ice Cap retreated by 31 m while the southeastern margins advanced by up to 29 m (Bol'shiyanov and Makeev, 1995). It is believed that the migration of the Vavilov Ice Cap in a southwesterly direction is due to a shift in the accumulation zone toward the present moisture source, following changes in the general atmospheric circulation over the region (Barkov et al., 1992; Bol'shiyanov and Makeev, 1995). This notion is supported by palynological analysis of an ice core extracted from the Vavilov Ice Cap (Andreev et al., 1997). The presence of considerable amounts of pollen from a west European species of lime in the upper 65 m of the core suggests that summer air masses started to be predominantly from the southwest approximately 500 years ago.

The aim of this paper is to comprehend the flow dynamics of the Vavilov Ice Cap in the context of the present climate and to test the hypothesis that the ice cap is not in equilibrium with the present climate. To do this, a mass balance model of the Vavilov Ice Cap was coupled to an ice-flow model and run under conditions representative of those at present (Bassford et al. 2006). The results quantify, for the first time, the flow dynamics of the ice cap and their relationship with the climate of this remote archipelago.

Mass Balance Model

The mass balance model used here has been adapted from the model of Greuell and Konzelmann (1994), and full details are available in Bassford (2002). The model calculates the surface energy balance, the temperature of the snow and ice cover, and rates of melting, refreezing, and run-off. Bassford et al. (2006) use the model to establish the mass balance of the Vavilov Ice Cap at present, using meteorological data from Vavilov Station as input (Fig. 2c). The refreezing of meltwater either directly on the ice surface (superimposed ice) or within the snowpack (refrozen ice) accounts for the bulk of the positive accumulation of ice under the present climate. Clearly, if refreezing processes were absent, the ice cap would be in significant negative balance.

The Ice-flow Model

CONTINUITY EQUATION

The basic equation at the heart of the ice-flow model is the continuity equation for ice (Mahaffy, 1976). This equation relates time dependent changes in ice thickness to the net ice flux and mass balance. With horizontal axes x and y, the vertically integrated continuity equation may be written as

where (ū.h) and (v̄.h) are the fluxes of ice in the x and y directions, respectively, and are calculated by the product of ice thickness, h, and the x or y component of the depth averaged velocity, ū or v̄. The mass balance b includes the processes of surface accumulation and ablation, basal melting and mass loss through iceberg calving. Since the latter two processes do not apply to the Vavilov Ice Cap, b is equivalent to the surface mass balance and is calculated by the distributed mass balance model (Bassford et al., 2006).

INTERNAL ICE DEFORMATION

Solving the full force balance for all of the stress components operating in a glacier is a complex problem. A common approach to simplify ice-flow modeling is the shallow ice approximation (Hutter, 1983). This assumes that ice deforms by shearing in horizontal planes parallel to the ice surface, known as laminar flow, and that shear stresses are the only non-zero stress component in the force balance. However, longitudinal stresses may become important in ice flowing over an undulating bed or in parts of a glacier with large spatial gradients in basal motion. Neither of these conditions applies to the Vavilov Ice Cap. Therefore, the shallow ice approximation is a valid approach.

For laminar flow, the relation between the shear strain rate ε̇_xz and the shear stress τ_xz can be expressed using the flow law proposed by Glen (1955):

where u is the x-component of velocity. In this study, the terms A and n, the flow law parameter and exponent, are equal to 8.47 × 10⁻¹⁶ s⁻¹ (kPa)⁻³ and 3, respectively. The shear stress at depth (h − z) is obtained by

Here ρ_i is the density of ice, g is acceleration due to gravity (9.81 m s⁻²), and θ is the surface gradient in the x direction (Paterson, 1994). For small values of θ, as is generally the case for the surface gradients of an ice cap, sin θ may be approximated by θ. Substituting this simplified expression into Equation 3 and integrating with respect to z gives

where u_s is the surface velocity and u is the velocity at depth (h − z) due to internal deformation. Assuming that basal motion is negligible for the ice cap in this study, the depth-averaged velocity ū is

The same series of equations applies to calculate the velocity component in the y direction, v̄. Having derived the depth-averaged velocity in terms of ice cap thickness and surface elevation, ū and v̄ can be substituted back into Equation 1 to solve for time dependent changes in ice thickness.

The Coupled Mass Balance and Ice-flow Model

Interaction between the mass balance and ice-flow models is important as several feedback mechanisms have a powerful effect on the response of a glacier to climate change. For example, the elevation–mass balance feedback mechanism is the process through which changes in specific balance are amplified by changes in surface elevation. Thus, a long-term negative specific balance causes a lowering of the ice surface which in turn leads to increased ablation due to higher temperatures at a lower elevation. The feedback works in reverse for a positive mass balance. A different process which acts against the elevation–mass balance feedback is the interaction between ice cap hypsometry and mass balance. If a glacier has a long-term positive mass balance and is advancing, then the flux of ice will increase across the equilibrium line leading to a larger ablation area which reduces the mass balance. The opposite applies to a retreating glacier.

The majority of models used to examine the sensitivity of glacier mass balance to climate change assume that the glacier geometry remains fixed in time. This sensitivity is based on the present hypsometry of a glacier and is termed the static mass balance sensitivity to climate change (Warrick et al., 1996). As the geometry of a glacier adjusts to new mass balance conditions, the change in mass balance over the varying area and elevation is the dynamic sensitivity. The latter represents the complete solution to predicting how a glacier may respond to climate change. Although calculating the dynamic sensitivity of glaciers is difficult on a regional scale (Oerlemans et al., 1998), it is possible for individual ice masses and is the objective in this study.

Since the computation of the surface mass balance is quite time demanding, especially when calculations are made for each of the grid cells in the model domain, it is important to consider how frequently the mass balance should be redefined in the ice-flow model. At the first time step, mass balance is calculated for every cell in the grid, including cells outside the margins of the initial ice cap. Subsequent computing time can be reduced by avoiding repeat mass balance calculations that are unnecessary in parts of the ice cap which have changed very little since the last mass balance calculation. This is particularly relevant to the Vavilov Ice Cap, which evolves very slowly over time. Therefore, the models are coupled in such a way that the ice-flow model continues to run until the surface elevation of a grid cell changes by more than a certain amount since the last time at which mass balance was calculated for that cell. Mass balance is then redefined for the cell, before the ice-flow model resumes. This means that mass balance calculations are concentrated in the parts of the ice cap that are changing the most. The number of new mass balance calculations will decrease to zero as the coupled model reaches a steady state.

The mass balance model is called when there is a change in elevation of ±50 m; equating to a difference in mass balance of only about ±5 cm w.e. (water equivalent) a⁻¹. Uncertainties in other parts of the model do not warrant more frequent recalculations of the surface mass balance.

Present-day Flow Dynamics of the Vavilov Ice Cap

Figure 4 shows that the velocity field calculated by the model for the present day Vavilov Ice Cap is characterized by low surface velocities, generally less than 5 m a⁻¹. The distribution of ice velocity follows a classical flow structure for a dome-shaped ice cap, in which velocity increases from zero at the ice divide to a maximum close to the equilibrium line, before decreasing toward the ice margins (Paterson, 1994). An anomaly to this general pattern is the relatively fast flow in the southwestern part of the ice cap, where surface velocities increase to 9 m a⁻¹ (Fig. 4). The relatively fast flow at the margin is expected since this part of the ice cap is currently advancing (Fig. 3). In agreement with the model results, ice surface velocities of 0.5–4.5 m a⁻¹ have been measured in the northern part of the ice cap, with greater velocities of 5–25 m a⁻¹ recorded at the southern margin (Barkov et al., 1992).

A synthetic aperture radar interferogram, based on images acquired in spring 1996, shows the velocity field over the entire Vavilov Ice Cap in the west–east look direction of the ERS-1 satellite (Fig. 5). The absence of interference fringes between the ice divide and margins around the majority of the ice cap implies low surface velocities of generally <10 m a⁻¹ in the look direction of the satellite. However, an area of relatively fast flow in the western part of the ice cap is clearly evident from the more closely spaced fringes, which indicate that velocity increases to a maximum of about 20 m a⁻¹ at the ice margin (Fig. 5). The pattern of fringes shows that the fast flow unit extends back into the slow moving ice to within about 5 km of the ice divide.

The effect of the fast flow unit on the interior of the ice cap is manifested in the shape of the ice surface elevation contours which show a slight trough at the head of the fast flow unit. Similar fast flow features have been reported for other Arctic ice caps, including the Academy of Sciences Ice Cap in Severnaya Zemlya (Dowdeswell et al., 2002) and Austfonna in Svalbard (Dowdeswell et al., 1999).

Comparison of Figures 4 and 5 highlights two major differences between the modeled distribution of surface velocity over the Vavilov Ice Cap and the flow structure determined from SAR interferometry. First, the model calculates an area of relatively fast flow in the southwestern part of the ice cap which is not apparent from the interferogram. This disparity could occur because interferometry is not sensitive to surface displacement perpendicular to the look direction of the satellite. Therefore, the velocity of ice flowing toward the southwestern margin is poorly sensed by the satellite and not well represented in the interferogram. Second, the fast flow unit in the western part of the ice cap identified by interferometry is not simulated by the model. Interferometrically derived surface velocities of 3–20 m a⁻¹ for the fast flow unit are significantly greater than the modeled surface velocities of only 1–2 m a⁻¹. This suggests that ice motion due to internal deformation is augmented by a considerably larger basal component in the western part of the ice cap, assumed to be negligible in the model. Since much of the ice cap bed in the marginal zone of the fast flow unit lies close to or below sea level (Fig. 2b), it is possible that basal motion may involve the deformation of overridden marine sediments. It should be noted that the upper part of the fast flow unit is located in a subglacial valley, supporting the notion that bed topography is an important influence on the location of these ice stream-like features (Dowdeswell et al., 2002). In other parts of the ice cap, where both modeling and interferometry indicate slow surface velocities generally <5 m a⁻¹, a frozen bed condition is likely, and ice motion occurs solely through internal deformation.

Time Dependent Simulations

Two time dependent simulations were performed in which the model was run for 1000 years to examine the evolution of the Vavilov Ice Cap under the present climate. In the first simulation, referred to hereafter as the “reference run”, the original assumption was maintained that ice motion due to basal sliding is negligible for the entire ice cap. In a second simulation, referred to as the “fast flow run”, special treatment was given to the western side of the ice cap by forcing the model to reproduce the interferometrically derived surface velocities of the fast flow unit (Fig. 5). This was achieved implicitly by augmenting the velocity calculated from internal deformation by a basal motion component, so that the total surface velocity was equal to that determined from interferometry. Velocity in other parts of the ice cap was calculated from internal deformation alone, as in the reference run. The distribution of surface velocity for the present-day Vavilov Ice Cap calculated by this scheme is shown in Figure 6.

Time dependent changes in the geometry of the Vavilov Ice Cap for the reference run are shown in Figure 7. During the 1000 year period, the northern and western margins of the ice cap gradually retreat by an average of 5 km, the eastern and southern margins remain stable, while the southwestern part of the ice cap advances by up to 5 km. Overall, there is a reduction in ice cap area of 170 km², equivalent to a decrease of 10% in the original area, although the total volume of ice increases by 20 km³, about 3% of the initial volume (Fig. 8). The latter is reflected by a thickening of ice in the center of the ice cap and an increase in the surface elevation of the summit from 708 to 752 m a.s.l. (Fig. 7). The general migration of the ice cap toward the southwest is consistent with the observed changes in ice cap geometry. However, there is a noticeable difference between the modeled and measured positions of the western margin, which is currently advancing (Fig. 8).

As expected, incorporating the fast flow unit in the model has a large effect on the evolution of the western part of the Vavilov Ice Cap (Fig. 9). In the first 200 years of the fast flow run, the model calculates a relatively rapid advance of about 1.5 km in the western margin at the terminus of the fast flow unit, resulting in an initial increase of ∼10 km² in the total area of the ice cap (Fig. 8). An advance of the magnitude calculated by the model would cause this part of the ice cap to terminate in marine waters, although iceberg calving was not accounted for in the model. The fast flow unit has the effect of draining ice from the interior of the ice cap toward the margins, reflected by a dip in surface elevation extending back to within 5 km of the ice divide (Fig. 9). This leads to a greater initial reduction in ice cap volume relative to the reference run, as more ice is transported into the ablation zone (Fig. 8). After 200 years it appears that the fast flow unit becomes exhausted of ice, resulting in a retreat of the western margin as it ablates without the extra flux of mass from the accumulation zone. This implies that the relatively fast flow rates observed in the western part of the ice cap are not sustainable under the present climate. Throughout the 1000 year simulation, the fast flow unit has little impact on the behavior of other parts of the ice cap which evolve in a similar way as in the reference run.

In both the reference and fast flow runs, the Vavilov Ice Cap failed to reach a steady state, suggesting that (1) the ice cap is not in equilibrium with the present-day climate, and (2) the ice cap has a response time with an order of magnitude of 10³ years. It is, therefore, possible that the ice cap is still responding to changes in climate identified by ice core analysis to have occurred about 500 years ago (Andreev et al., 1997).

Summary and Conclusions

A coupled mass balance and ice-flow model was used to calculate the velocity field over the Vavilov Ice Cap. In agreement with ground based measurements, the model indicates that ice flow in the bulk of the ice cap occurs through internal deformation alone, producing low surface velocities generally <5 m a⁻¹. However, an area of relatively fast flow in the western part of the ice cap, where surface velocities reach 20 m a⁻¹, was identified from SAR interferometry. It is likely that motion in the fast flow unit is increased by a significant basal component. Time dependent simulations suggest that the fast flow unit is responsible for the observed advance of a section of the western margin of the ice cap, although it appears that the relatively fast flow rates measured by interferometry are not sustainable for more than about 200 years under the present climate. Model results are consistent with recent observations which show that the ice cap is currently migrating toward the precipitation source in the southwest, supporting the hypothesis that the ice cap is not in equilibrium with the present climate. The response time of the Vavilov Ice Cap is estimated to have an order of magnitude of 10³ years, making it possible that the ice cap is still responding to changes in climate thought to have occurred about 500 years ago.